# Modelo das camas no hospital em Banks et al, 2012

a <- c(
	"U*(m*nu1+beta*(C+J*(1-gamma)))",
	"J*(m*nu2)",
	"C*((1-m)*nu2)",
	"J*((1-m)*nu2)",
	"C*(alpha)")

V <- t(matrix(c(
	-1,+1, 0,
	 0,+1,-1,
	+1,-1, 0,
	+1, 0,-1,
	 0,-1,+1), ncol=3, byrow=TRUE))

parms <- c(
	m = 0.04,      # infected admission rate
	beta = 0.001,  # effective contact rate
	gamma = 0.58,  # hand hygiene compliance rate
	alpha = 0.29,  # patient isolation rate
	nu1 = 0.16,    # uninfected discharged rate
	nu2 = 0.08     # infected discharged rate
)

# X( #Uncolonized, #Colonized, #J-isolation )
N <- 37  #beds
x0 <- N*c(U=29/37,C=4/37,J=4/37)

library(GillespieSSA)
library(gridExtra)
library(ggplot2)
library(reshape2)

theme_set(theme_gray(base_size=12) %+replace% theme(
	axis.text        = element_text(color='black'),
	panel.grid.minor = element_line(color='gray90'),
	panel.grid.major = element_line(color='gray60'),
	panel.background = element_rect(fill='white'),
	legend.text.align=0,
	legend.key       = element_blank(),
	legend.position  = 'bottom'
))

g <- list()
for(method in c('D','ETL','BTL','OTL')) {
	out <- ssa(x0, a, V, parms, 200, method)
	colnames(out$data)[1] <- 'time'
	rownames(out$data) <- 1:nrow(out$data)
	df <- as.data.frame(out$data)
	df <- melt(df, id='time', value.name='beds', variable.name='type')
	g <- c(g, list(ggplot(df, aes(x=time, y=beds, color=type)) + geom_line() + ggtitle(method) + scale_x_continuous("time (days)")))
}

# http://stackoverflow.com/questions/17150183/r-plot-multiple-lines-in-one-graph

p <- do.call(arrangeGrob, g)
ggsave('plots.pdf', p)

df <- data.frame()
for(method in c('ETL','BTL','OTL'))  # c('D','ETL','BTL','OTL')
	for(N in 37*10^(0:2)) {
		print(N)
		x0 <- N*c(U=29/37,C=4/37,J=4/37)
		t <- system.time(ssa(x0, a, V, parms, 200, method))[1]
		df <- rbind(df, data.frame(N=N, method=method, time=t))
	}

p <- ggplot(df, aes(x=N, y=time, color=method)) + geom_line() + geom_point() + scale_x_log10() + scale_y_continuous("time (seconds)") + theme(axis.title.x = element_blank())
ggsave('times2.pdf', p, width=8, height=3)
